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ABSTRACT 

Recent experiments [Nakata, M. ef a/., End-to-end 
stacking and liquid crystal condensation of 6 to 20 
basepair DNA duplexes. Science 2007; 318:1276- 
1279] have demonstrated spontaneous end-to-end 
association of short duplex DNA fragments 
into long rod-like structures. By means of exten- 
sive all-atom molecular dynamic simulations, 
we characterized end-to-end interactions of duplex 
DNA, quantitatively describing the forces, free 
energy and kinetics of the end-to-end association 
process. We found short DNA duplexes to spontan- 
eously aggregate end-to-end when axially aligned 
in a small volume of monovalent electrolyte. It 
was observed that electrostatic repulsion of 
5 -phosphoryl groups promoted the formation of 
aggregates in a conformation similar to the B-form 
DNA double helix. Application of an external force 
revealed that rupture of the end-to-end assembly 
occurs by the shearing of the terminal base pairs. 
The standard binding free energy and the kinetic 
rates of end-to-end association and dissociation 
processes were estimated using two complemen- 
tary methods: umbrella sampling simulations of 
two DNA fragments and direct observation of the 
aggregation process in a system containing 458 
DNA fragments. We found the end-to-end force to 
be short range, attractive, hydrophobic and only 
weakly dependent on the ion concentration. The 
relation between the stacking free energy and 
end-to-end attraction is discussed as well as 
possible roles of the end-to-end interaction in 
biological and nanotechnological systems. 

INTRODUCTION 

Self-assembly properties of nucleic acids are vital to the 
basic functions of a biological cell and have been exten- 
sively exploited in biotechnology. DNA hybridization — 



self-assembly of complementary sequence single-stranded 
DNA (ssDNA) into a double helix — is a central biotech- 
nological process (1), used, among others, in platforms for 
DNA detection (2), programmable assembly of DNA 
nanostructures (3,4), directional transport of cargo (5), 
molecular computing (6) and nanofabrication (7,8). 
Another process of outstanding importance is DNA con- 
densation, where counterions transform electrostatic 
repulsion between naked DNA molecules into attraction, 
facilitating packaging of double-stranded DNA (dsDNA) 
in cell nuclei and viral capsids (9,10). 

Recently, an entirely different type of DNA self- 
assembly was discovered: spontaneous end-to-end aggre- 
gation of short duplex DNA fragments into rod-like struc- 
tures (11). When water was evaporated from solution 
containing a high concentration of short (6-20 bp) DNA 
fragments, liquid crystal phases were observed. Since the 
DNA fragments were nearly as wide as they were long, the 
observation of axial ordering could only be explained if 
the fragments formed rod-like supramolecules, suggesting 
end-to-end aggregation. Further experimental evidence of 
end-to-end association was obtained from the analysis of 
small angle X-ray scattering data from a system contain- 
ing short DNA fragments and a divalent electrolyte 
(12,13). The second virial coefficient extracted from 
these data was shown to be positive for DNA fragments 
capped with a short hairpin (indicating overall repulsion) 
and negative for DNA fragments without such caps 
(indicating overall attraction). It was concluded that 
end-to-end attraction was large enough to overcome elec- 
trostatic repulsion in a divalent electrolyte. While 
side-by-side forces between long DNA molecules has 
been the subject of many experimental (14-16) and theor- 
etical (17,18) studies, little is known about the conditions 
and microscopic mechanism of DNA association end-to- 
end. Furthermore, the effects of end-to-end attraction of 
duplex DNA in biological and technological processes are 
entirely unexplored. 

Whereas traditional single molecule experiments have 
provided extensive information about DNA hybridization 
and side-by-side interactions (14,19), applying these tools 
to study end-to-end assembly is extremely difficult, as a 
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DNA duplex cross section is just 5 nm 2 and the effective 
concentration of DNA ends in a solution amenable to 
single-molecule studies is typically small. All-atom 
molecular dynamics (MD) method is well suited for the 
study of systems that share the length scale of short 
dsDNA molecules and can be used to probe the atomic 
origin of intermolecular forces (20). Here, we use the MD 
method to characterize end-to-end association of duplex 
DNA in unprecedented detail, elucidating the microscopic 
mechanism of spontaneous association, its free energy 
costs and the kinetic rates. At the end of the article, we 
discuss the relationship between our findings and the per- 
tinent experimental observations. 

MATERIALS AND METHODS 

General simulation methods 

All MD simulations were performed using the program 
NAMD (21), the parmbscO refinement of the AMBER- 
parm99 force field (22,23), the TIP3P model of water 
(24), standard parameters for ions (25), periodic 
boundary conditions, and particle-mesh Ewald (PME) 
full electrostatics with a PME grid density of ~1A per 
grid point. Except where specified, van der Waals and 
short-range electrostatic energies were calculated using a 
smooth (10-12 A) cutoff, and integration was performed 
using 1-2-4 fs multiple timestepping (21). The temperature 
was held constant using a Langevin thermostat (21) applied 
to all non-hydrogen atoms; the Langevin damping constant 
was set to 0.1 /ps. For simulations in the NPT ensemble, 
constant pressure was maintained at 1 bar using the 
Nose-Hoover Langevin piston pressure control (26). 

Each simulation reported in this study used one of the 
following three system types: elongated along the z-axis to 
minimize the amount of solvent around two DNA frag- 
ments (^24000 atoms, Figure la); isotropic to allow two 
DNA fragments to tumble freely (~56 000 atoms, Figure 
2a); and large and isotropic to allow unbiased interaction 
between 458 DNA fragments (~1.4 M atoms, Figure 4a). 
The DNA sequence was poly(dA-dT) in all systems. 
Counterions were added to each system to neutralize the 
DNA charge prior to the addition of a number of ions 
corresponding to the reported molarity (100 mM, except 
where specified) of NaCl electrolyte. Steric clashes that 
were introduced during the assembly of each system 
were removed from each system through minimization 
using a conjugate gradient method (27). Equilibration 
was performed in the NPT ensemble, and subsequent pro- 
duction simulations were performed in the NVT ensemble, 
except where specified. 

Collapse of aligned dsDNA 

Thirty-six systems were built using the anisotropic unit 
cell. The dimensions of each system were chosen to 
provide a minimum distance of 2nm between the 
surfaces of the DNA fragments across the periodic 
boundary, which should accommodate a majority of 
screening counterions for NaCl solutions of lOOmM or 
greater concentration (Debye length <lnm). Steric 
clashes were removed through 3000 minimization steps. 



Each system was subsequently equilibrated for 65 ps 
with the DNA backbone atoms harmonically restrained 
to their initial positions. Axial alignment of the DNA frag- 
ments was enforced by harmonically restraining each 
phosphorous atom of the DNA to the surface of an 
11 -A radius cylinder (with spring constants of 139 pN/ 
nm per atom). The DNA fragments could translate 
along the axis of the cylinder and rotate azimuthally^. 
The starting conformation was characterized by a 20.5 A 
end-to-end separation, which we define as the distance 
between the centers of mass of the nearest terminal base 
pairs, taking the periodic boundary condition into 
account. The relative azimuthal angle § of the terminal 
base pairs was defined as the angle between the projections 
of the vectors connecting the 05' and 03' atoms of the 
terminal base pairs into the plane normal to the common 
DNA axis. For two consecutive base pairs in a B-DNA 
helix, 4> « 36°, depending on the sequence. 

Stability of the end-to-end complex 

Systems were built by placing collapsed end-to-end DNA 
assemblies in an isotropic volume of lOOmM NaCl elec- 
trolyte (Figure 2a). Steric clashes with the solvent were 
removed through minimization having the DNA 
backbone restrained. Subsequent simulation was per- 
formed in the NPT ensemble. 

Mechanics of end-to-end dissociation 

The coordinates of the system containing 
S'-phosphorylated DNA fragments in a lOOmM electro- 
lyte were taken after 140 ns of simulation described in the 
section 'Stability of the end-to-end complex' to provide 
the initial conformation for simulations of rupture of the 
DNA fragments. The NVT ensemble was used during 
these simulations; Langevin thermostat was applied to 
water oxygens and ions. A harmonic spring of £ s = 4000 
pN/nm was used to produce the rupture, by increasing its 
rest length at a rate of 0.4 or 0.2 A/ns. The work per- 
formed at both rates were in good agreement. 
Simulations of the end-to-end fragments having 
complementary single-stranded overhangs are described 
in Section 1.4 of Supplementary Data. 

Potential of mean force of axially aligned DNA duplexes 

Umbrella sampling simulations were performed using the 
anisotropic systems (Figure la) and two simulation proto- 
cols different by the method used to set up initial systems 
and the alignment restraints. Both protocols enforced the 
end-to-end distance r using a harmonic spring of 
£ s = 4000 pN/nm for 3.5 < r < 12 A o in 0.5-A intervals 
and k s = 1000 pN/nm for 13 < r < 19 A in 1.0-A intervals. 
The first protocol was used to provide the estimate of the 
potential of mean force (PMF) (Figure 3). The initial con- 
formation for each simulation was obtained by placing the 
DNA fragments a specified distance r apart at one of the 
four 4> = 0, 90, 180 and 270° (four simulations for each r). 
The systems were equilibrated for 2 ns before data accu- 
mulation during production simulations lasting M6ns. In 
the second protocol, which we used to compute the 
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relative binding free energies (Table 1), the initial con- 
formations were generated iteratively by shifting the 
minimum of the restraining potential in steps and 
followed by 0.5-ns equilibration, starting from the final 
frames obtained in the simulations described in the 
section 'Collapse of aligned dsDNA'. Subsequently, each 
system was equilibrated for at least 2.5 ns before accumu- 
lation of data during production simulations lasting 
7.5-1 5 ns. Axial alignment was maintained as described 
in the section "Collapse of aligned dsDNA", using 
k s = 13.9 and 139 pN/nm for the first and second proto- 
cols, respectively. In the second protocol, a torque 
pointing along the common DNA axis was distributed 
among the phosphorous atoms of each DNA molecule 
to restrain <\> about —20°, 36° or 180°, with a spring 
constant of 219.4 pNnm/rad 2 , which roughly corresponds 
to an 8° root mean squared fluctuation. 

Spontaneous assembly of long end-to-end aggregates 

The system (depicted in Figure 4a) was assembled through 
the sequential placement of 458 DNA fragments in a cubic 
volume (250 A on each side) of lOOmM NaCl electrolyte. 
To place a DNA fragment, trial positions and orientations 
were randomly selected until the DNA coordinates did not 
clash with any previously placed fragments. During the 
first 50 ns of equilibration, the system shrank to its equi- 
librium size of 238 A on each side. To improve computa- 
tional efficiency, a 7-8 A cutoff was used along with 2-2-6 
fs time stepping scheme. 

RESULTS 

Collapse of aligned dsDNA 

Spontaneous end-to-end association of duplex DNA was 
observed in the simulations of two (dA-dT) 10 fragments 



constrained to diffuse along a common axis in a volume of 
lOOmM NaCl electrolyte. Figure la illustrates the initial 
state of a typical simulation system. Figure lb plots the 
distance between the DNA fragments versus time for two 
simulation systems differed by the termination of the 
DNA's Spends. The DNA fragments were observed to 
freely diffuse along the common axis until the 
end-to-end distance fell below « 8 A, whereupon the frag- 
ments rapidly collapsed into an end-to-end bound 
complex. The relative azimuthal orientation of the DNA 
fragments § continued to change after the collapse, 
switching between several metastable orientations 
(Figure lc). We define the relative azimuthal angle c|) as 
the angle between the projections of the vectors connect- 
ing the 05' and 03' atoms of the terminal base pairs into 
the plane normal to the common DNA axis (see 'Materials 
and Methods' section). 

In the final conformation adopted by the blunt-ended 
fragments, the 5' to 3' direction of the backbone was dis- 
continuous at the end-to-end junction. In the case of the 
S'-phosphorylated fragments the 5' to 3' direction of the 
backbone was continuous at the end-to-end junction as 
though the backbone of a continuous 20 bp B-DNA had 
been cut. 

To determine the statistical significance of the above 
observations, we performed 17 additional simulations 
for each system type, different by the relative azimuthal 
orientation of the DNA fragments at the onset of the 
simulation: <fy t = 0 = ix20 o , where /=1,..., 17. In all 
simulations, we observed collapse of the DNA fragments 
into an end-to-end bound complex. Figure Id plots the 
relative azimuthal orientation at the time of collapse 
against the time to collapse. The collapse of blunt-ended 
fragments occurs irrespective of their azimuthal orienta- 
tion, whereas formation of a S'-phosphorylated end-to- 
end assembly did not occur with 4> between 90 and 230°. 




time (ns) time (ns) 

Figure 1. Collapse of aligned dsDNA. (a) Simulation system containing axially aligned duplex DNA. Each DNA fragment is free to move along and 
rotate about the common axis. The DNA duplexes (blue and green) are shown in van der Waals representation; sodium and chloride ions are shown 
as yellow and cyan spheres; water is not depicted. An animation illustrating spontaneous end-to-end collapse of duplex DNA is available in 
Supplementary Data, (b and c) The end-to-end distance (b) and the relative azimuthal angle § (c) of two duplex DNA in representative simulations 
of the end-to-end collapse. Data from the same pair of simulations are plotted in (b) and (c). (d and e) Scatter plot showing the relative azimuthal 
angle <$> at the time of collapse (d), and at the end of simulation (e). One data point is shown for each of 36 simulations of blunt-ended (black circles) 
and 5 ; -phosphorylated (red squares) dsDNA fragments in 100 mM NaCl electrolyte. 
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Figure 2. Stability of the end-to-end DNA assembly, (a) Simulation system containing two spatially unrestrained dsDNA fragments in 100 mM 
NaCl electrolyte. Both fragments are free to rotate and move about the simulation box. The system is drawn as in Figure la. An animation 
illustrating a typical MD trajectory is available in Supplementary Data, (b and c) End-to-end junction of S'-phosphorylated dsDNA with § = —20° 
(b) and blunt-ended dsDNA with 4> = 180° (c). (d and e) End-to-end distance (d) and relative azimuthal angle, (e) of the DNA fragments in three 
unrestrained MD trajectories. 



After collapse, 4> continued to evolve, reaching the 
states characterized in Figure le. The clustering of 4> 
values around —20, 36 and 180° indicates the three 
preferred binding states. At <\> « 36°, the conformation 
of the bound complex is similar to that of a continuous 
B-form DNA. At $ « -20°, the backbones of the 
terminal base pairs overlap slightly (Figure 2b). At 
4> ~ 180°, the 5'-ends of the fragments are in direct 
contact (Figure 2c). The conformations of the 
end-to-end assemblies with 4> = —20 and 36° are similar, 
differing primarily in the order in which the 5'- and 
3 / -termini are lapped and in their base stacking 
geometry. Thus, relative to the coordinates of two con- 
secutive base pairs in a canonical B-DNA helix, the 
terminal base pairs forming an end-to-end junction have 
a time-averaged root mean squared deviation (RMSD) of 
2.1 and 0.9 A for the § = —20 and 36° conformations, 
respectively. For reference, base pairs in the middle of 
one of the DNA fragments had an RMSD of 0.7 A. 

The preference for these three orientations suggests a 
hydrophobic origin of the attractive force, as such con- 
formations reduce exposure of the hydrophobic bases to 
water. The relative orientations of the blunt-ended DNA 
were nearly equally split between the —20 and 180° states. 
More than 50% of the S'-phosphorylated fragments 
formed the —20° state, 35% formed either the 36° or 
180° state (three systems each) and 12% formed the 
state with § « 100°. We attribute such preferential align- 
ment of the S'-phosphorylated fragments to the electro- 
static repulsion between the terminal phosphate groups, 
which apparently renders the 180° orientation energetic- 
ally less favorable than the —20° one. The free energy 
difference between these bound states is discussed below. 



Stability of the end-to-end complex 

DNA fragments initially forming a bound state were 
simulated in the absence of any restraints using the iso- 
tropic system shown in Figure 2a. Three systems were 
constructed: one containing S'-phosphorylated DNA 
bound with <fy = — 20° (Figure 2b) and two containing 
blunt-ended DNA fragments bound at $ = 180° 
(Figure 2c) and 4> = —20°. All three systems contained 
lOOmM NaCl electrolyte. 

The plot of the end-to-end distance, Figure 2d, indicates 
that all three assemblies remained bound for the entire 
duration of the simulations (> 200 ns). The standard devi- 
ation of the end-to-end distance in the S'-phosphorylated 
system was 0.68 A, twice less than that of the blunt-ended 
systems. The greater stability of the S'-phosphorylated 
complex may be due to hydrogen bonds that were 
observed between the S'-terminal phosphate of one 
DNA fragment and the S'-terminal hydroxyl of the other 
fragment. The plot of the relative azimuthal angle reveals 
that the S'-phosphorylated system maintained the 
same stable conformation at § « —20°, depicted in 
Figure 2b. Starting from a similar conformation, the 
blunt-ended complex underwent two sudden rotations at 
70 and 140 ns that brought <|> from -20° to -110° and to 
180°; the relative orientation continued to evolve after 
that. 

In the simulation of the 180° blunt-ended complex, the 
terminal base pair of one of the fragments ruptured after 
33 ns. During the next few nanoseconds, Watson-Crick 
pairs within that fragment stochastically broke and 
re-annealed, propagating the unpaired base toward the 
opposite end of the DNA fragment, and slipping the 
entire DNA strand with respect to the other by 1 bp. 
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Despite the slippage, the DNA fragments remained stably 
bound. 

Mechanics of end-to-end dissociation 

The lifetime of a bound complex sharply decreases when 
an external force is applied to disrupt it (28). Thus, under 
a constant force of 150pN directed along the common 
axis of the DNA fragments, the assembly remained 
intact during a 50-ns simulation but dissociated within 
5 ns under a 200 pN force (see Section 1.1 in 
Supplementary Data for simulation details). 

To determine the dissociation pathway, the rupture 
force and the mechanical work required to dissociate the 
assemblies, the DNA fragments were subjected to the 
force of a harmonic spring whose equilibrium-extension 
length was increased at a constant rate. Spring-driven 
rupture of this sort has been used extensively in the 
study of proteins (29,30). Rupture was induced by 
pulling apart the fragments either along (axial stretching) 
or perpendicular to (transverse shearing) their common 
symmetry axis, applying the force either to the nearest 
ends or to the centers of mass (CoM) of the fragments. 
At least four simulations were performed for each 
protocol (14 in total). Animations in Supplementary 
Data illustrate typical simulation trajectories. In all 
cases, rupture was observed to occur by sliding of one 
terminal base pair relative to the other and was preceded 
by stretching of the duplex in the case of the CoM axial 
pulling. Although the three rupture protocols yielded dif- 
ferent dependencies of the force on the separation distance 
(Supplementary Figures SI and S2), the average work 
performed was 9.4 ± 1.5 kcal/mol, independent of the 
rupture protocol. The typical rupture forces varied 
between 100 and 200 pN and were considerably lower 
for CoM pulling. Inclusion of short overhangs of compli- 
mentary sequence ssDNA at the ends of the fragments 
increased the work required to rupture the end-to-end 
assemblies (for details, see Section 1.4 and Figure S3 in 
Supplementary Data). 

PMF of axially aligned DNA duplexes 

To improve our estimates of the force and free energy of 
the end-to-end interaction, 100 variants of the system 
containing two-DNA fragments in lOOmM NaCl 
electrolyte, Figure la, were simulated with a constant 
end-to-end distance enforced by a harmonic spring poten- 
tial. Each simulation explored a unique combination of 
the end-to-end distance and the azimuthal angle and 
lasted 18 ns (aggregate simulation time was 1.86 jis). The 
DNA fragments were kept aligned by weak harmonic re- 
straints (see 'Materials and Methods' section) that allowed 
the terminal base pairs to shear. 

Figure 3 shows the dependence of the effective 
end-to-end force on the end-to-end distance and the 
PMF reconstructed from this set of simulations by the 
weighted histogram analysis method (31,32). The PMF 
can be thought of as the change of free energy along a 
chosen coordinate. The force sharply increases with the 
end-to-end distance between 3.5 A — the distance between 
consecutive base pairs in a DNA helix — and 6.5 A, the 



separation allowing water molecules to penetrate the 
volume between the ends of the fragments. The force 
rapidly decreases as the end-to-end distance exceeds 
6.5 A and becomes slightly repulsive after ~13A (2-10 
pN). Thus, the end-to-end force has a large but very 
short-range attractive component caused by the hydro- 
phobic effect and a much smaller long-range repulsive 
component that originates from screened electrostatic 
interactions between the DNA fragments (19). 

A variation of the above protocol (described in Section 
2 of Supplementary Data) was used to calculate the 
PMF for DNA fragments different by their terminal 
chemistry and relative azimuthal orientation and for 
several concentrations of the surrounding electrolyte. 
Table 1 lists the change in the depth of the PMF 
minima (A(min[PMF] - PMF(oo)) relative to its value 




J i l i l i l 

5 10 15 20 

end-to-end distance (A) 



Figure 3. Representative dependence of the effective force (red) and the 
free energy (black) of two axially aligned DNA fragments on the 
end-to-end distance. The data result from 100 independent simulations 
of two S'-phosphorylated fragments, the end-to-end distance of which 
was maintained at a specified value by a harmonic spring. Additional 
restraints were applied to maintain the axial alignment. The strength of 
the restraints was found to affect the values but not the general shape 
of both curves (see text). The image in the background illustrates the 
simulation method. The DNA fragments were immersed in 100 mM 
NaCl electrolyte, and were free to rotate about their helical axes. 



Table 1. Relative free energy change, A AG, upon formation of the 
end-to-end complex 



5' chemistry 


ion concentration 
(M) 




A AG (kcal/mol) 


Phosphoryl 


0.1 


36° 


0 


Phosphoryl 


1.0 


36° 


0.4 


Phosphoryl 


1.0 


-20° 


-1.8 


Phosphoryl 


1.0 


180° 


0.1 


Hydroxyl 


1.0 


36° 


2.3 


Hydroxyl 


1.0 


180° 


1.2 



Here, the free energy change AG is approximated by the minimum of 
the end-to-end PMF obtained from umbrella sampling simulations 
(Supplementary Figure S4). The values of AG are given relative to 
the value measured for the system containing S'-phosphorylated frag- 
ments end-joined in the 36° orientation in 100 mM NaCl. In each 
simulation, the relative azimuthal orientation of the DNA fragments 
was enforced using harmonic restraints (see 'Material and Methods' 
section). The application of such restraints introduced a bias to the 
estimates of AG. As all simulations employed the same restraints, 
this bias was assumed to cancel out in the calculation of AAG. 
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for the S'-phosphorylated fragments with § = 36° in 
lOOmM NaCl. These calculations (detailed in Supple- 
mentary Data) demonstrate that increasing the electrolyte 
concentration from 0.1 to 1 M has negligible effect on the 
PMF. Among the three most likely orientations that 
S'-phosphorylated fragments form, the § = —20° state 
has the lowest free energy, in good agreement with the 
4> = —22° angle frequently observed in the crystal struc- 
tures of poly(dA-dT) oligonucleotides (33). The S'-phos- 
phorylated fragments exhibit deeper minima than the 
blunt-ended fragments. For the latter, the conformation 
of 4> = 180° is preferred over 4> = 36°. These variations in 
the depth of the PMF are consistent with the occupancy of 
bound states observed in our simulations of spontaneous 
collapse of aligned DNA fragments (Figure le). 

Standard binding free energy of end-to-end association 

We computed the standard free energy of binding G bind 
using a variation of the method described previously (34). 
For a system of two DNA fragments, G b i nd determines the 
fraction of time the fragments form a bound state, which 
can be, in principle, observed directly in an all-atom MD 
simulation. However, the dissociation rate of the 
end-to-end assembly was found to exceed 200 ns 
(Figure 2d) and therefore using such a brute force 
approach was not possible. Equivalently, G bind can be 
obtained from the logarithm of the equilibrium binding 
constant, which is the ratio of the kinetic rates of end 
joining and rupture (k on and fc off , respectively). 

For our calculations of G bind , we considered a process 
consisting of the following four steps (Supplementary 
Figure S5). First, we evaluated the free energy cost of 
enforcing axial alignment restraints on a pair of infinitely 
separated DNA fragments. Second, we computed the cost 
of bringing a pair of DNA fragments from infinity to the 
maximum-separation state considered in our simulations 
(CoM-CoM distance of 52 A). Third, we determined the 
free energy cost of forming the end-to-end complex of two 
axially aligned DNA fragments. Finally, we evaluated the 
cost of releasing the axial alignment restraints from 
the end-to-end assembly. The sum of these terms yielded 
the free energy change upon formation of the end-to-end 
DNA complex G bind = — 6.3 ± 1 kcal/mol for a DNA con- 
centration of 1 M in 120mM NaCl. Since a pair of DNA 
ends has only one binding configuration, we can express 
the standard binding free energy in terms of DNA ends, so 
that G^nd = —5.4 kcal/mol. A complete description of the 
methods used in this section is provided in Section 3.1 in 
Supplementary Data. 

The above value for G b md represents our best effort to 
quantify the strength of the end-to-end interaction. In 
general, G bind may depend on the DNA sequence. 
Although the question of sequence dependence is 
enticing, we defer investigations of the sequence depend- 
ence of the end-to-end interaction to future studies. 

The rate of dissociation of an end-to-end assembly, £ of f, 
can be estimated from the PMF and diffusion coefficient, 
D, by computing the mean first passage time under the 
assumption that reannealing does not occur (35,36) (see 
Section 3.2 in Supplementary Data for details). From the 



simulations performed in the section 'Collapse of aligned 
dsDNA', we estimate D ~25 A 2 /ns. Due to uncertainty in 
the location of the barrier peak, we obtain the range 
k~l -170 000-860 000 ns, with a best estimate of 
k~ff —480 000 ns. The use of axial alignment restraints in 
our calculations of the PMF limits pathways ordinarily 
available to rupturing DNA, so these values likely repre- 
sent an upper bound for k~#. 

Assuming that end-to-end binding is limited by transla- 
tional diffusion, the upper bound estimate for k on is 
4kDR 0 ^1 (ns M) _1 , where R 0 = 37 A is the CoM 
distance between a pair of DNA fragments at which 
binding is expected to occur. The true value of k on 
should be smaller because the DNA fragments must be 
axially aligned for binding to occur. Furthermore, 
long-range electrostatic repulsion may reduce the value 
of k on . Our estimates of G b md and fc off suggest a range 
for k on of 0.03-0.16 ns" 1 M. 

Spontaneous assembly of long end-to-end aggregates 

Our results until this point described the interaction of two 
DNA fragments in isolation. To simulate multifragment 
aggregation, 458 DNA fragments, each 10 bp in length, 
were placed in a cube of lOOmM NaCl solution (23.8 nm 
on each side) to form the system shown in Figure 4a. 
During a 260-ns MD simulation, the DNA fragments 
diffused about their initial positions and interacted with 
their neighbors to form aggregates up to 11 DNA frag- 
ments (1 10 bp) in length. The DNA fragments that formed 
the longest 10 aggregates are shown at the beginning, 
Figure 4b, and the end, Figure 4c, of the simulation. 
The number of aggregates of a given length at three dif- 
ferent instances of the MD trajectory is shown in 
Figure 4d. The plot reveals rapid growth of end-to-end 
aggregates and roughly exponential distribution of the 
lengths of the aggregates. The number of aggregates did 
not reach a steady state by the end of the simulation 
(Supplementary Figure S6). 

We model the process of end-to-end aggregation using a 
simplified reversible step-growth polymerization model, 
which has an analytical time-dependent solution that 
relates the mean aggregation number, (AO to G bind as a 
function of DNA concentration, c (37) (see Section 4 in 
Supplementary Data). The DNA concentration in our 
multifragment system was ~56mM. Using our prediction 
of Gbind = — 6.3 kcal/mol, we obtain (TV) = 39 for the 
mean aggregation number at equilibrium for this system. 

The above simulation partially mimics the experimental 
assay of Nakata et al. (11), which involved very dense 
fluids of short DNA fragments. For 10 bp DNA frag- 
ments, Nakata et al. found a transition from an isotropic 
phase to a nematic liquid crystal phase at c = 875 mg/ml 
(~150mM). Through a combination of calibration experi- 
ments and theory, the authors estimated (TV) = 9 at the 
isotropic-nematic phase transition. Under the model of 
reversible step growth polymerization, (AO =9 at 
c = 150mM corresponds to G bind = —3.87 kcal/mol. In 
contrast, using our prediction of G bind = —6.3 kcal/mol 
yields (AO = 64 at a DNA concentration of 150 mM. 



3818 Nucleic Acids Research, 2012, Vol. 40, No. 9 




(d) 



100=- 



10| 



1=- 



MD: 4 ns 
MD: 32 ns 
MD: 256 ns 

RSGP: G = -3.9 kcal/mol 
RSGP: G = -6.3 kcal/mol 



4 8 
length of aggregate, s 



12 



Figure 4. Spontaneous aggregation of duplex DNA into long rod-like 
structures, (a) A system containing 458 duplex DNA fragments placed 
at random. NaCl solution is shown as a semi-transparent molecular 
surface, (b and c) Initial (b) and final (c) conformations of the DNA 
fragments that composed the longest 10 aggregates at the end of a 
260-ns MD simulation. DNA fragments forming each aggregate are 
shown in a different color, (d) The instantaneous number of aggregates, 
N(s), formed by s DNA fragments in the MD simulation. The lines 
show the equilibrium distribution of the aggregate length according to 
the reverse step-growth polymerization model (37) for specified values 
of Gbind- A movie of the simulation trajectory is provided in 
Supplementary Data. 

In the reversible step growth polymerization model, the 
total number of associated molecules is described by a 
second-order kinetic equation. Since unbinding of DNA 
fragments is negligible in our system, k on can be extracted 
by fitting N(t) = No /(I + tc 0 k on ) to the data, where N(t) is 
the total number of end-to-end bound DNA fragments at 
time t, N 0 = N(0) and c 0 is the initial DNA concentration. 
In our simulation of spontaneous aggregation, the associ- 
ation rate reduced from k on = 0.37 ns/M observed within 
the first 50 ns to 0.069 ns/M for the 60-250-ns interval 
(Supplementary Figure S6). The change in association 
rate occurred as longer end-to-end aggregates formed 
(>3 duplex fragments) and the rotational and translation- 
al diffusive motions of shorter aggregates slowed. 

The aggregation simulation provides a lower bound 
estimate for the dissociation rate k Q ff. During the course 
of the simulation, 307 DNA ends were bound for 217 ns 
on average and unambiguous unbinding was observed for 
only one end-to-end associated complex. A few events 
were observed in which partial unbinding occurred via 



rupture of the Watson-Crick base pair of terminal nucleo- 
tides, as well as one instance where a bound DNA 
fragment was transferred to an unbound fragment; we 
neglect these events for the subsequent analysis. A 
Poisson distribution yields the probability that exactly 
one end-to-end bound DNA pair would rupture during 
the simulation given a value of k Q ff. For k~^< 1900 ns, it 
is likely that more than one rupture would have occurred. 
For k~ff > 59600 ns, it is likely that no rupture would have 
occurred. Thus, the probability that exactly 1 rupture 
occurred is >10% for = 19 000-596 000 ns. The 
greatest probability of exactly 1 rupture occurring is 
37% for k~ff = 67 000 ns. The ranges estimated for k on 
and k 0 ff suggest — Gbind m between 4.4 and 7.6 kcal/mol, 
with a most probable range of 5.2-6.2 kcal/mol 
[k on = 0.069-0.37 ns/M and k~^ = 67, 000ns]. Alternative 
statistical analysis of the results yields a similar range for 
G bind and is described in Section 4 of Supplementary Data. 



DISCUSSION AND CONCLUSIONS 

Stability of a dsDNA molecule can be conveniently 
described as a sum of base stacking and base pairing 
energies contributed by individual nucleotides that 
constitute the molecule. It is tempting to conceptually 
equate the base stacking interactions within a continuous 
molecule with the base stacking interactions that drive the 
assembly of two disjoint DNA duplexes. Thus, the unified 
nearest neighbor parameters, which can predict the energy 
for DNA hybridization based on DNA melting data, 
suggest G bin d = -16.94 + 2 x 6.94 = -3.06 kcal/mol (38) 
for the association of two 10 bp poly(dA-dT) DNA frag- 
ments into a continuous 20-bp molecule. Such a simple 
calculation may be, however, flawed as additional con- 
formational flexibility afforded by the lack of phospho- 
diester bonds at the end-to-end interface should allow 
the base stacking geometry to be optimized, magnifying 
the base stacking contribution to the free energy. 
Accordingly, the average interaction energy between 
adjacent base pairs with § = —20° was measured to be 
2 kcal/mol lower in bases forming the end-to-end 
junction than in bases in the middle of one of the DNA 
fragments (Supplementary Figure S7). This finding is in 
good agreement with a survey of crystallographic struc- 
tures that found DNA fragments formed of AT base pairs 
to stack with 4> = —22° (33). We note that the stacking 
geometry may depend on the sequence of the DNA 
termini. 

The base stacking free energy has been experimentally 
quantified by introducing a dangling nucleotide or a nick 
(a cut in the backbone of one strand) to a DNA molecule 
and observing the change in melting temperature (39-41), 
and by observing the mobility of a nicked DNA molecule 
relative to intact DNA and DNA with a gap (42,43). These 
experiments provide estimates for the base stacking free 
energy between —0.65 and —2.0 kcal/mol for stacks 
formed by thymine and adenine. However, the extraction 
of the free energy values is indirect with these methods. 

Here, we provide the first direct estimate of the standard 
binding free energy of end-to-end association of DNA 
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fragments, obtaining G b md = —6.3 ± 1 kcal/mol. The 
simulations reported here are similar to the liquid crystal 
condensation experiments of the Clark and Bellini groups 
(11), which estimated —3.8 kcal/mol for the end-to-end 
free energy, in reasonable agreement with our estimate. 
The small-angle X-ray scattering experiments of the 
Pollack group demonstrated that the end-to-end inter- 
action dominates over electrostatic repulsion in a 
divalent electrolyte (12,13), which, unfortunately, is not 
sufficient to estimate the standard free energy of 
end-to-end binding. We note that the value of Gbind 
obtained in this study is larger than values reported in 
experiments. Having employed multiple methods, we are 
confident that the range of values obtained for Gbind ac- 
curately reflects the standard binding free energy within 
the limits of the molecular force field used in our study. 
Nevertheless, we cannot rule out the possibility that the 
present MD force field somewhat exaggerates the inter- 
actions driving end-to-end self-assembly of duplex DNA. 

It is interesting to note that hydrophobic interactions 
between DNA bases and inorganic materials can be sig- 
nificantly stronger than the stacking energies observed in 
biochemical assays. AFM experiments indicate that 
ssDNA adheres to graphite with free energies of —4.9 
and —6.8 kcal/mol per nucleotide for cytosine and 
guanine, respectively (44). For comparison, a computa- 
tional investigation of hydrophobic interactions revealed 
a free energy of —55 kcal/mol for the adhesion of a pair of 
11 x 12 A 2 graphene sheets, or about 10 times the energy 
per unit surface area (45). Interactions of similar strength 
were found to promote DNA-fullerene and ssDNA- 
carbon nanotube association (46,47). 

Given the relatively large free energy of the end-to-end 
interaction, we pose the following question: why has 
end-to-end association only recently been observed? In 
most biological and nanotechnological settings, the concen- 
tration of DNA ends is too low for end-to-end association 
to be statistically significant, and the lifetime of end-to-end 
interactions (~70 |is) is somewhat shorter than the temporal 
resolution of many experimental techniques. To illustrate 
this point, we plot in Figure 5 the fraction of bound DNA 
ends versus the concentration of free DNA ends in chemical 
equilibrium. The concentration of DNA ends is relatively 
high in DNA cyclization assays that measure the fraction of 
DNA molecules bent into a circle. The /-factor, which rep- 
resents the concentration of one end in the proximity of the 
other, has a maximum value of M0 -4 mM (48), which is 
two orders of magnitude too small to promote cyclization of 
a significant fraction of blunt-ended DNA. The introduc- 
tion of sticky ends increases the interaction energy by 
a few kcal/mol so that the DNA cyclization can be observed. 

The large standard binding energy of end-to-end asso- 
ciation implies that end-to-end interactions will be import- 
ant in systems containing a high concentration of DNA 
ends. For instance, a remarkable new method of creating 
patterned, self-assembled structures out of DNA — termed 
DNA origami (3) — introduces many nicks along a path of 
DNA, which may enhance the stability of the resultant 
pattern. We speculate that broadened awareness of 
end-to-end association will influence the development of 
nanotechnologies where the end-to-end interaction can be 




of DNA ends (M) 



Figure 5. The effect of end-to-end attraction in different DNA systems. 
The binding free energy (black) and the fraction of bound DNA ends 
(red) are plotted against the reference concentration of DNA ends. 
Background images schematically illustrate four DNA systems in 
which the end-to-end attraction may or may not play a role. From 
top left to bottom right: the maximum local concentration of DNA 
ends [i.e. the /-factor (51)] is about two orders of magnitude too low to 
induce an observable fraction of blunt-ended DNA circles of any length 
(orange); translation and rotational confinement of DNA ends at 
dsDNA breakage [e.g. by the protein Ku (PDB:1JEY)] will promote 
binding of the DNA ends, which likely aids repair of DNA during 
non-homologous end joining (49) (purple); the structure factor 
obtained from small-angle X-ray scattering experiments of short 
DNA duplexes in a divalent electrolyte reveals end-to-end attraction 
(12) (blue); at very high DNA concentrations, long DNA aggregates 
form and align in liquid crystal phases (11) (green). 

used advantageously, such as with DNA origami, or can 
pose a limitation, such as for DNA microarrays. 

In cells, double-stranded DNA breakage, which poses a 
mortal threat, results in nearby blunt or sticky DNA ends. 
During non-homologous end-joining — the dominant 
repair pathway for such breaks in multicellular eukaryotes 
(49,50) — the ruptured DNA ends are held together by 
proteins such as the Ku heterodimer or DNA-PK, or by 
nucleosome interactions until damaged DNA can be 
removed and ligation of the DNA backbone occurs (49). 
Since DNA attracts end-to-end, it is not necessary that 
this complex, whose microscopic structure is not yet 
known, hold the DNA ends in strict alignment. It is suf- 
ficient for the ends to be held proximally so that the effect- 
ive concentration of DNA ends is large enough to 
promote end-to-end association (Figure 5, purple), where- 
upon ligation may occur. The free energy of dsDNA 
end-to-end association found in this study suggests that 
placing DNA ends in a sphere of 3-nm radius will produce 
an end-to-end associated state ~95% of the time. 

SUPPLEMENTARY DATA 

Supplementary Data are available at NAR Online: 
Supplementary Methods, Results and Figures 1-7, 
Supplementary Animations 1-6 and Supplementary 
References [52—61]. 
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